home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / External Modules / External modules sources / Pascal / GuideLine.p < prev    next >
Text File  |  1996-04-26  |  24KB  |  743 lines
  1. {******************************************************************************}
  2. { GuideLine.p                                                            }
  3. {                                                                         }
  4. {                                                                         }
  5. { Version 22.04.96                                                    }
  6. {******************************************************************************}
  7.  
  8.  
  9. unit GuideLine;
  10.  
  11. interface
  12. {$IFC UNDEFINED THINK_PASCAL }
  13.     uses
  14.         Types, Memory, proFit_interface;
  15. {$ELSEC}
  16.     uses
  17.         proFit_interface;
  18. {$ENDC}
  19.  
  20.     procedure main (selector: integer; pb: ExtModulesParamBlockPtr);
  21.  
  22.  
  23.  
  24. implementation
  25.  
  26. {$IFC NOT UNDEFINED OLDROUTINENAMES }
  27.     {$IFC NOT OLDROUTINENAMES }
  28.     procedure DisposPtr (p: Ptr);    { if we don't have defined DisposPtr }
  29.     begin
  30.         DisposePtr(p);
  31.     end;
  32.     {$ENDC}
  33. {$ENDC}
  34.  
  35. { note: MPW users must make sure that the procedure main is at the beginning of the compiled code }
  36. { under Think Pascal, this is cared for by the compiler }
  37. { We let main call a function mainMain to make sure that the code starts with a jump to }
  38. { our entry point even when compiling under MPW Pascal }
  39.     procedure mainMain (selector: integer; pb: ExtModulesParamBlockPtr);
  40.     forward;
  41.     procedure main (selector: integer; pb: ExtModulesParamBlockPtr);
  42.  
  43.     begin
  44.         mainMain(selector, pb);
  45.     end;
  46.  
  47.  
  48.  
  49.     const
  50.         maxGuidePts = 100;        { maximum number of points to define the smooth curve }
  51.  
  52.     type
  53.         MyGlobals = record
  54.                 nrGuidePts: integer;        { effective number of points defining the smooth curve }
  55.                 xv, yv, mv: array[1..maxGuidePts] of extended;{ points and slope of this curve }
  56.             end;
  57.         ptrToMyGlobals = ^MyGlobals;
  58.  
  59.         ValsT = array[1..8192] of real;
  60.         ValsP = ^ValsT;
  61.  
  62.  
  63.     function CubicCurve (xm, x1, y1, m1, x2, y2, m2: extended): extended;
  64. { 2 points with their slope given define a cubic polygon; the y value at xm is calculated. }
  65.         var
  66.             aa, bb, cc, dd: extended;
  67.             dx, dy, x1x1: extended;
  68.     begin
  69.         dx := x1 - x2;
  70.         dy := y1 - y2;
  71.         aa := (dx * (m1 + m2) - 2 * dy) / (dx * dx * dx);
  72.         x1x1 := x1 * x1;
  73.         bb := ((m1 - m2) - 3 * aa * (x1x1 - x2 * x2)) / 2 / dx;
  74.         cc := m1 - 3 * aa * x1x1 - 2 * bb * x1;
  75.         dd := y1 - aa * x1x1 * x1 - bb * x1x1 - cc * x1;
  76.         CubicCurve := aa * xm * xm * xm + bb * xm * xm + cc * xm + dd;
  77.     end;        { CubicCurve }
  78.  
  79.  
  80.     procedure FindAveragedVal (nrPoints: integer; xvals, yvals: ValsP; log: boolean;    {}
  81.                                     AA, reg, eps, xxx: extended; var yyy: extended);
  82. { this function calculates an averaged value at x=xxx, by summing over all points given, }
  83. { and weigthing them linearly or logarithmically and using either linear or quadratic regression. }
  84.         var
  85.             s1, s2, s3, s4, s5, s6, s7, s8: extended;    { summing parameters }
  86.             xi, yi, gi: extended;                            { coordinates and weight of point i }
  87.             m, b: extended;                                    { parameters for linear regression }
  88.             va, vb, vc: extended;                            { parameters for quadratic regression }
  89.             val1, val2: extended;                            { found for averaging }
  90.             ex, ex1: extended;
  91.             i, ccount: integer;
  92.     begin
  93.         s1 := 0;        { sum gi   (weights of points xi) }
  94.         s2 := 0;        { sum gi*xi }
  95.         s3 := 0;        { sum gi*yi }
  96.         s4 := 0;        { sum gi*xi*yi }
  97.         s5 := 0;        { sum gi*xi*xi }
  98.         s6 := 0;        { sum gi*xi*xi*xi }
  99.         s7 := 0;        { sum gi*xi*xi*xi*xi }
  100.         s8 := 0;        { sum gi*xi*xi*yi }
  101.         ccount := 0;
  102.  
  103.         for i := 1 to nrPoints do
  104.             begin
  105.                 if TestStop then
  106.                     exit(FindAveragedVal);
  107.                 xi := xvals^[i];
  108.  
  109.                 if log then                            { different weighting factors for log and lin cases }
  110.                     gi := exp(-sqr(AA * ln(xi / xxx)))
  111.                 else
  112.                     gi := exp(-sqr(AA * (xi - xxx)));
  113.  
  114.                 if gi > 1e-9 then
  115.                     begin
  116.                         yi := yvals^[i];
  117.                         ccount := ccount + 1;
  118.  
  119.                         if ccount <= 3 then                { if two points have the same x-value ... }
  120.                             if ccount = 1 then
  121.                                 val1 := xi
  122.                             else if ccount = 2 then
  123.                                 if abs(xi - val1) < eps then
  124.                                     ccount := 1
  125.                                 else
  126.                                     val2 := xi
  127.                             else { if ccount = 3 then }
  128.                                 if (abs(xi - val1) < eps) | (abs(xi - val2) < eps) then
  129.                                     ccount := 2;
  130.  
  131.                         s1 := gi + s1;                        { summing }
  132.                         s2 := gi * xi + s2;
  133.                         s3 := gi * yi + s3;
  134.                         s4 := gi * xi * yi + s4;
  135.                         ex := gi * xi * xi;
  136.                         s5 := ex + s5;
  137.                         if reg = 2 then
  138.                             begin
  139.                                 s6 := ex * xi + s6;
  140.                                 s7 := ex * xi * xi + s7;
  141.                                 s8 := ex * yi + s8;
  142.                             end;
  143.                     end;
  144.             end;    { for i := 1 to nrPoints do }
  145.  
  146. { final evaluation: }
  147.         if ccount = 0 then
  148.         else if ccount = 1 then                        { no regression }
  149.             yyy := s3 / s1
  150.         else
  151.             begin
  152.                 ex := s5 * s1 - s2 * s2;
  153.                 if (reg = 1) | (ccount = 2) then            { linear regression }
  154.                     begin
  155.                         m := (s4 * s1 - s2 * s3) / ex;
  156.                         b := (s3 - m * s2) / s1;
  157.                         yyy := m * xxx + b;
  158.                     end
  159.                 else
  160.                     begin                                    { quadratic regression }
  161.                         ex1 := s5 * s2 - s6 * s1;
  162.                         va := (ex * (s8 * s1 - s3 * s5) - ex1 * (s3 * s2 - s4 * s1)) / (ex * (s7 * s1 - s5 * s5) - sqr(ex1));
  163.                         vb := (va * ex1 + (s4 * s1 - s3 * s2)) / ex;
  164.                         vc := (-va * s5 - vb * s2 + s3) / s1;
  165.                         yyy := va * xxx * xxx + vb * xxx + vc;
  166.                     end;
  167.             end; { if count >1 then }
  168.     end;    { FindAveragedVal }
  169.  
  170.  
  171.  
  172. {***********************************************************************************************}
  173.  
  174.     procedure SetUp (var moduleKind: integer;    { set moduleKind to isFunction or isProgram }
  175.                                     var name: Str255;             { the name of the program or function }
  176.                                     var requiredGlobals: longint;     { the number of bytes to be allocated in ExtModulesParamBlock.globals }
  177.                                         { set requiredGlobals to 0 if you don't use this feature }
  178.                                     pb: ExtModulesParamBlockPtr);    { the complete parameter block passed by pro Fit to the }
  179.                                         { routines defined in this file. In most cases it can be ignored }
  180. { SetUp is called once when the external module is linked to pro Fit }
  181.     begin
  182.         moduleKind := isFunction;
  183.         name := 'GuideLine';
  184.         requiredGlobals := sizeof(MyGlobals);
  185.     end;
  186.  
  187.  
  188. {***********************************************************************************************}
  189.  
  190.     procedure InitializeFunc (var hasDerivatives: boolean;    { set this to true if and only if you define the function }
  191.                                                 { Derivatives to calculate the partial derivatives of the parameters }
  192.                                     var descr1stLine, descr2ndLine: Str255;     { The two lines of the text in the parameter window }
  193.                                     var numberOfParams: integer;             { the number of parameters of the function }
  194.                                     var a0: DefaultParamInfo;                 { the default names, values etc. of the parameters }
  195.                                     pb: ExtModulesParamBlockPtr);            { the complete parameter block passed by pro Fit to the }
  196.                                                                     { routines defined in this file. In most cases it can be ignored }
  197. { InitializeFunc is called once (after SetUp has been called) when the external module is linked to pro Fit }
  198. { Used to set all the information needed to describe a function }
  199.     begin
  200.         hasDerivatives := false;
  201.  
  202.         descr1stLine := 'smoothness: 0-100 %; interval: 0=auto, 1=linear, 2=logarithmic, 3=point density; ';
  203.         descr2ndLine := '$BSelect data$averaging: 0=mean, 1=linear regr., 2=quadratic regr.';
  204.                         { Using $B...$ in descr2ndLine creates a button in the parameter window, }
  205.                         { which, when pressing, calls the function Check() with 30000 as argument. }
  206.  
  207.         numberOfParams := 3;
  208.         a0.value^[1] := 80;
  209.         a0.mode^[1] := constant;
  210.         a0.lowest^[1] := 0;
  211.         a0.highest^[1] := 100;
  212.         a0.name^[1] := 'smooth';
  213.  
  214.         a0.value^[2] := 0;
  215.         a0.mode^[2] := constant;
  216.         a0.lowest^[2] := 0;
  217.         a0.highest^[2] := 3;
  218.         a0.name^[2] := 'interval';
  219.  
  220.         a0.value^[3] := 0;
  221.         a0.mode^[3] := constant;
  222.         a0.lowest^[3] := 0;
  223.         a0.highest^[3] := 2;
  224.         a0.name^[3] := 'average';
  225.  
  226.         pb^.v[1] := 0.31938153;
  227.         pb^.v[2] := -0.356563782;
  228.         pb^.v[3] := 1.781477937;
  229.         pb^.v[4] := -1.821255978;
  230.         pb^.v[5] := 1.330274429;
  231.         pb^.v[6] := 0;        { the x-column to be used, 0 if not defined yet. }
  232.         pb^.v[7] := 0;        { the y-column to be used, 0 if not defined yet. }
  233.         pb^.v[8] := 0;        { the window ID of the selected data window, 0 if not defined yet. }
  234.     end;
  235.  
  236.  
  237. {***********************************************************************************************}
  238.  
  239.     function Check (ParamNo: integer;             { the parameter that was changed }
  240.                                     var a0: DefaultParamInfo;         { the default names, values etc. of the parameters }
  241.                                     pb: ExtModulesParamBlockPtr    { the complete parameter block passed by pro Fit to the}
  242.                                         { routines defined in this file. In most cases it can be ignored }
  243.                                     ): CheckPAnswer;
  244.     { Can be left emtpy (returning good) if not needed. }
  245.     { called when the user has changed a value in the parameters window. This routine }
  246.     { can then check if this parameters is fine. It can also change some of the }
  247.     { other entries in a0. The returned values can be: }
  248.     {    good:        return this value if you agree with the new parameter value }
  249.     {    update:        return this value if you want the parameters window }
  250.     {                to be updated because you changed some of the values in a0 }
  251.     {    bad:        return this value if you want the new parameter value to be refused }
  252.         var
  253.             inprec: inputRec;                                { variables to call InputBox() }
  254.             inpID, inpXcol, inpYcol, pHelp: extended;        { data window id and column number of x- and y-values }
  255.             inpIDStr, inpXcolStr, inpYcolStr, pHelpStr: Str255;    { texts for inputBox }
  256.     begin
  257.         Check := good;
  258.         if (ParamNo = 2) & (a0.value^[2] <> 0) & (a0.value^[2] <> 1) & (a0.value^[2] <> 2) & (a0.value^[2] <> 3) then
  259.             Check := bad;
  260.         if (ParamNo = 3) & (a0.value^[3] <> 0) & (a0.value^[3] <> 1) & (a0.value^[3] <> 2) then
  261.             Check := bad;
  262.         if (ParamNo = 30000) then
  263.             begin
  264.                 if (pb^.v[6] = 0) then
  265.                     inpXcol := xColumn        { generates an error if there is no data window }
  266.                 else
  267.                     inpXcol := pb^.v[6];
  268.                 if (pb^.v[7] = 0) then
  269.                     inpYcol := yColumn
  270.                 else
  271.                     inpYcol := pb^.v[7];
  272.                 if (pb^.v[8] = 0) then
  273.                     inpID := GetCurrentWindow(dataType)
  274.                 else
  275.                     inpID := pb^.v[8];
  276.  
  277.                 pHelp := 0;
  278.                 inpIDStr := '$WData window :';
  279.                 inpXcolStr := '$CX column :';
  280.                 inpYcolStr := '$CY column :';
  281.                 pHelpStr := '$XPrint help into Results window';
  282.                 inprec[1].x := @inpID;
  283.                 inprec[1].s := @inpIDStr;
  284.                 inprec[2].x := @inpXcol;
  285.                 inprec[2].s := @inpXcolStr;
  286.                 inprec[3].x := @inpYcol;
  287.                 inprec[3].s := @inpYcolStr;
  288.                 inprec[4].x := @pHelp;
  289.                 inprec[4].s := @pHelpStr;
  290.                 SetBoxTitle('GuideLine Setup');
  291.                 if InputBox(4, inprec) then
  292.                     begin
  293.                         StopExecution;
  294.                         exit(Check);
  295.                     end;
  296.                 pb^.v[6] := inpXcol;
  297.                 pb^.v[7] := inpYcol;
  298.                 pb^.v[8] := inpID;
  299.                 if pHelp > 0 then
  300.                     begin
  301.                         WritelnString('');
  302.                         WritelnString('');
  303.                         WritelnString('Short description of the function GuideLine:');
  304.                         WritelnString('-------------------------------------------');
  305.                         WritelnString('GuideLine generates a curve going smoothly through a given data set.');
  306.                         WritelnString('');
  307.                         WritelnString('The curve is found in the following way:');
  308.                         WritelnString('a) The given data set is divided into subsets of data points.');
  309.                         WritelnString('b) From each subset a "synthetic data point" is generated.');
  310.                         WritelnString('   The synthetic data point is the average of all points in');
  311.                         WritelnString('   the subset.');
  312.                         WritelnString('c) GuideLine returns a cubic spline curve going through');
  313.                         WritelnString('   the synthetic data points.');
  314.                         WritelnString('');
  315.                         WritelnString('The number of subsets is adjusted by the parameter smooth:');
  316.                         WritelnString(' - a[1] = 0  : there will be 100 subsets.');
  317.                         WritelnString(' - a[1] = 100: there will 1 subset only.');
  318.                         WritelnString('The algorithm for selecting the subsets (intervals) is');
  319.                         WritelnString('defined by a[2]:');
  320.                         WritelnString(' - a[2] = 0: automatic.');
  321.                         WritelnString(' - a[2] = 1: linear (subsets have same width on linear scale).');
  322.                         WritelnString(' - a[2] = 2: logarithmic (subsets have same width on log. scale).');
  323.                         WritelnString(' - a[2] = 3: point density (subsets have same number of data).');
  324.                         WritelnString('The averaging method is defined by a[3]: ');
  325.                         WritelnString(' - a[3] = 0: mean values of x- and y-coordinates.');
  326.                         WritelnString(' - a[3] = 1: y = linear regression at the mean x-coordinate.');
  327.                         WritelnString(' - a[3] = 2: y = quadratic regression at the mean x-coordinate.');
  328.                         WritelnString('');
  329.                     end;
  330.                 Check := good;
  331.             end;
  332.     end;        { Check }
  333.  
  334.  
  335. {***********************************************************************************************}
  336.  
  337.     procedure First (a: ParamArray;             { the new parameters }
  338.                                     pb: ExtModulesParamBlockPtr);    { the complete parameter block passed by pro Fit to the}
  339.                                         { routines defined in this file. In most cases it can be ignored }
  340.  
  341. { Can be left emtpy if not needed. }
  342. { Called whenever the parameters were changed. Can be used to accelerate }
  343. { some calculations. See manual for more info }
  344.         label
  345.             10;
  346.         var
  347.             f2max, f2min, f2mean: extended;    { extremal values of given data points }
  348.             interval: extended;                    { interval width, within which a smooth point is calculated }
  349.             eps: extended;                        { smallest distance between the x-values of two data points }
  350.             lowB, upB: extended;                { boundaries of interval as x-values }
  351.             lowN, upN: integer;                { boundaries of interval as index of data points }
  352.             xvals, yvals: ValsP;                { arrays containing the given data points }
  353.             AA: extended;                        { exp. weighting factor }
  354.             xmean, xd: extended;                { center of weighting }
  355.             ymean, yd: extended;                { y vals }
  356.             nrInts: extended;                    { number of intervals }
  357.             nrPoints: integer;                    { number of points given }
  358.             nrGuide: integer;                    { number of points calculated for smooth curve }
  359.             regr: integer;                        { regression: 0=mean value, 1=linear, 2=quadratic }
  360.             xcol, ycol: integer;                { the column numbers of the current data window }
  361.             i, j, k, count: integer;                { counting variables }
  362.             log: boolean;                        { true if 'logarithmic' weighting is used }
  363.             numint: boolean;                    { true if interval depends on point density }
  364.             P: ptrToMyGlobals;                { it is more comfortable to use P }
  365.             windID: longint;                        { window ID of data window to be used }
  366.     begin
  367.  
  368. {******************* prepare arrays ********************}
  369.         if (pb^.globals = nil) then
  370.             begin
  371.                 StopExecution;
  372.                 exit(First);
  373.             end;
  374.         P := ptrToMyGlobals(pb^.globals);
  375.  
  376.         if (pb^.v[6] = 0) then
  377.             xcol := xColumn        { generates an error if there is no data window }
  378.         else
  379.             xcol := round(pb^.v[6]);
  380.         if (pb^.v[7] = 0) then
  381.             ycol := yColumn
  382.         else
  383.             ycol := round(pb^.v[7]);
  384.         if (pb^.v[8] <> 0) then
  385.             begin
  386.                 windID := round(pb^.v[8]);
  387.                 SetCurrentWindow(windID);
  388.             end;
  389.  
  390.         numint := (a[2] = 3);
  391.         regr := round(a[3]);
  392.         xvals := ValsP(NewPtr(sizeof(ValsT)));
  393.         yvals := ValsP(NewPtr(sizeof(ValsT)));
  394.         if (xvals = nil) | (yvals = nil) then
  395.             begin
  396.                 if xvals <> nil then
  397.                     DisposPtr(Ptr(xvals));
  398.                 if yvals <> nil then
  399.                     DisposPtr(Ptr(yvals));
  400.                 if AlertBox('Not enough memory.') then
  401.                     ;
  402.                 StopExecution;
  403.                 exit(First);
  404.             end;
  405.         if TestStop then
  406.             exit(First);
  407.  
  408. {******************* get the points and sort if necessary ********************}
  409.  
  410.         f2mean := 0;
  411.         nrPoints := 0;
  412.         if (xcol = 0) | (ycol = 0) then
  413.             begin
  414.                 StopExecution;
  415.                 exit(First);
  416.             end;
  417.  
  418.         for i := 1 to NrRows do
  419.             if TestData(i, xcol) & TestData(i, ycol) then
  420.                 begin
  421.                     xd := GetData(i, xcol);
  422.                     j := nrPoints + 1;
  423.                     if numint & (nrPoints > 0) then        { sorting, if intervals dependent on point density }
  424.                         for j := nrPoints downto 1 do
  425.                             if xd > xvals^[j] then
  426.                                 begin
  427.                                     for k := nrPoints downto j + 1 do
  428.                                         begin
  429.                                             xvals^[k + 1] := xvals^[k];
  430.                                             yvals^[k + 1] := yvals^[k];
  431.                                         end;
  432.                                     j := j + 1;
  433.                                     goto 10;
  434.                                 end;        { for i }
  435. 10:
  436.                     xvals^[j] := xd;
  437.                     yvals^[j] := GetData(i, ycol);
  438.                     nrPoints := nrPoints + 1;
  439.                     if nrPoints = 1 then
  440.                         begin
  441.                             f2min := xd;
  442.                             f2max := xd;
  443.                         end;
  444.                     if f2min > xd then
  445.                         f2min := xd;
  446.                     if f2min > 0 then
  447.                         f2mean := f2mean + xd;
  448.                     if f2max < xd then
  449.                         f2max := xd;
  450.                 end;
  451.  
  452. {******************* check wether enough points ********************}
  453.  
  454.         if nrPoints < 1 then                                { at least 1 point }
  455.             begin
  456.                 StopExecution;
  457.                 if AlertBox('Not enough data for GuideLine.') then
  458.                     ;
  459.                 P^.nrGuidePts := nrPoints;
  460.                 StopExecution;
  461.                 exit(first);
  462.             end;
  463.         if nrPoints <= 2 then                            { easy cases }
  464.             begin
  465.                 P^.nrGuidePts := nrPoints;
  466.                 exit(first);
  467.             end;
  468.         if TestStop then
  469.             exit(First);
  470.  
  471. {******************* prepare variables ********************}
  472.  
  473.         nrInts := sqr((100 - a[1]) / 100) * (maxGuidePts - 1) + 1;
  474.                                     { a[1] is limited between 0 and 100 }
  475.                                     { nrInts will be between 1 and maxGuidePts }
  476.         nrGuide := 0;
  477.         eps := (f2max - f2min) / 2000;        { minimum to define two x-values to be different }
  478.         log := false;
  479.         if numint then
  480.             nrInts := round(nrInts);
  481.         interval := (f2max - f2min) / nrInts;
  482.         lowB := f2min - interval / 2;
  483.         if not numint & (f2min > 0) then
  484.             begin
  485.                 f2mean := f2mean / nrPoints;
  486. { logarithmic scale enforced or expected: }
  487.                 log := (a[2] = 2) | ((a[2] = 0) & (f2mean < (f2max - f2min) / 3.6));
  488.                 if log then
  489.                     begin
  490.                         interval := ln(f2max / f2min) / nrInts;
  491.                         lowB := exp(ln(f2min) - interval / 2);
  492.                     end;
  493.             end;
  494.         AA := 2 / interval;
  495.  
  496. {******************* calculate 1 point per interval ********************}
  497.  
  498.         i := 1;
  499.         while (numint & (i <= round(nrInts))) | (not numint & (lowB < f2max)) do
  500.             begin
  501.                 if numint then                        { calculate boundaries of interval }
  502.                     begin
  503.                         lowN := round((i - 1) / nrInts * nrPoints) + 1;
  504.                         upN := round(i / nrInts * nrPoints);
  505.                     end
  506.                 else if log then
  507.                     upB := exp(ln(lowB) + interval)
  508.                 else
  509.                     upB := lowB + interval;
  510.  
  511.                 xmean := 0;                            { sum up within interval }
  512.                 ymean := 0;
  513.                 count := 0;
  514.                 if numint then
  515.                     begin
  516.                         for j := lowN to upN do
  517.                             begin
  518.                                 xmean := xmean + xvals^[j];
  519.                                 ymean := ymean + yvals^[j];
  520.                                 count := count + 1;
  521.                             end;        { for j }
  522.                         i := i + 1;
  523.                     end
  524.                 else
  525.                     begin
  526.                         for j := 1 to nrPoints do
  527.                             begin
  528.                                 xd := xvals^[j];
  529.                                 if (xd >= lowB) & (xd < upB) then
  530.                                     begin
  531.                                         xmean := xmean + xd;
  532.                                         ymean := ymean + yvals^[j];
  533.                                         count := count + 1;
  534.                                     end;
  535.                             end;
  536.                         lowB := upB;
  537.                     end;
  538.  
  539.                 if count = 0 then                            { calculate mean value }
  540.                     cycle;
  541.                 nrGuide := nrGuide + 1;
  542.                 xmean := xmean / count;
  543.                 ymean := ymean / count;
  544.  
  545. {**************** calculate weight-averaged value and derivative ******************}
  546.  
  547.                 if regr > 0 then                                                    { calculate point 1 }
  548.                     FindAveragedVal(nrPoints, xvals, yvals, log, AA, regr, eps, xmean, ymean);
  549.                 P^.xv[nrGuide] := xmean;
  550.                 P^.yv[nrGuide] := ymean;
  551.                 if TestStop then
  552.                     exit(First);
  553.                 if regr > 0 then                                                    { calculate point 2 }
  554.                     begin
  555.                         if log then
  556.                             xd := exp(ln(xmean) + interval / 7)                        { delta x for slope ! (7 by accident) }
  557.                         else
  558.                             xd := xmean + interval / 7;
  559.                         yd := ymean;
  560.                         FindAveragedVal(nrPoints, xvals, yvals, log, AA, regr, eps, xd, yd);
  561.                         P^.mv[nrGuide] := (yd - ymean) / (xd - xmean);            { calculate the derivative }
  562.                     end;
  563.                 if TestStop then
  564.                     exit(First);
  565.             end;            { while (numint & i <= round(nrInts)) | (lowB < f2max) do }
  566.  
  567. {******************* clean up arrays********************}
  568.  
  569.         P^.nrGuidePts := nrGuide;
  570.         DisposPtr(Ptr(xvals));
  571.         DisposPtr(Ptr(yvals));
  572.  
  573. {******************* calculate slope for regr=0********************}
  574.  
  575.         if regr = 0 then
  576.             if nrGuide = 1 then
  577.                 P^.mv[1] := 0
  578.             else if nrGuide = 2 then
  579.                 begin
  580.                     P^.mv[1] := (P^.yv[2] - P^.yv[1]) / (P^.xv[2] - P^.xv[1]);
  581.                     P^.mv[2] := P^.mv[1];
  582.                 end
  583.             else
  584.                 begin
  585.                     for i := 2 to nrGuide do
  586.                         P^.mv[i] := (P^.yv[i] - P^.yv[i - 1]) / (P^.xv[i] - P^.xv[i - 1]);
  587.                     P^.mv[1] := P^.mv[2];
  588.                     for i := 2 to (nrGuide - 1) do
  589.                         P^.mv[i] := (P^.mv[i + 1] + P^.mv[i]) / 2;
  590.                     P^.mv[1] := 2 * P^.mv[1] - P^.mv[2];
  591.                     P^.mv[nrGuide] := 2 * P^.mv[nrGuide] - P^.mv[nrGuide - 1];
  592.                 end;
  593.     end;        {  First }
  594.  
  595.  
  596.  
  597. {***********************************************************************************************}
  598.  
  599.     procedure Func (x: extended;                 { the x-value }
  600.                                     a: ParamArray;                { the parameters }
  601.                                     var y: extended;             { the y-value }
  602.                                     pb: ExtModulesParamBlockPtr);    { the complete parameter block passed by pro Fit to the}
  603.                                         { routines defined in this file. In most cases it can be ignored }
  604. { called to calculate the y-value of your function for a given x and a given }
  605. { set of parameters }
  606.         var
  607.             P: ptrToMyGlobals;                { it is more comfortable to use P }
  608.             m, b: extended;                    { parameters for linear regression }
  609.             i, nrPoints: integer;
  610.     begin
  611.         if TestStop then
  612.             exit(Func);
  613.         if (pb^.globals = nil) then
  614.             begin
  615.                 StopExecution;
  616.                 exit(Func);
  617.             end;
  618.         P := ptrToMyGlobals(pb^.globals);
  619.  
  620.         nrPoints := P^.nrGuidePts;
  621.  
  622.         if nrPoints <= 1 then                { for the easy cases ... }
  623.             begin
  624.                 if nrPoints = 1 then
  625.                     y := P^.yv[1]
  626.                 else
  627.                     y := 0;
  628.                 exit(Func);
  629.             end;
  630.  
  631.         for i := 1 to nrPoints do            { find the two points, defining the curve piece needed }
  632.             if x < P^.xv[i] then
  633.                 leave;
  634.  
  635.         if i = 1 then                            { special case: curve left from the first point }
  636.             begin
  637.                 m := P^.mv[1];
  638.                 b := P^.yv[1] - m * P^.xv[1];
  639.                 y := m * x + b;
  640.             end
  641.         else if i = (nrPoints + 1) then        { special case: curve right to the last point }
  642.             begin
  643.                 m := P^.mv[nrPoints];
  644.                 b := P^.yv[nrPoints] - m * P^.xv[nrPoints];
  645.                 y := m * x + b;
  646.             end
  647.         else
  648.             begin                                { intermediate points }
  649.                 y := CubicCurve(x, P^.xv[i - 1], P^.yv[i - 1], P^.mv[i - 1], P^.xv[i], P^.yv[i], P^.mv[i])
  650.             end;
  651.     end;        { func }
  652.  
  653.  
  654.  
  655. {***********************************************************************************************}
  656.  
  657.     procedure Derivatives (x: extended;             { the x-value }
  658.                                     a: ParamArray;                 { the parameters }
  659.                                     var dyda: ParamArray;         { the derivatives }
  660.                                     pb: ExtModulesParamBlockPtr);    { the complete parameter block passed by pro Fit to the }
  661.                                         { routines defined in this file. In most cases it can be ignored }
  662.  
  663.     { Can be left empty if InitializeFunc sets hasDerivatives to false }
  664.     { called to calculate the partial derivatives of the function with respect to }
  665.     { its parameters. If you leave this function empty and set hasDerivatives to false in }
  666.     { FuncInitialize, the derivatives will be calcuated numerically, otherwise pro Fit }
  667.     { calls this function to obtain the values of ALL derivatives. }
  668.     { As a result of the numerical calculation fitting will be slower }
  669.     begin
  670.     end; { Derivatives }
  671.  
  672.  
  673.  
  674. {***********************************************************************************************}
  675.  
  676.     procedure Last (pb: ExtModulesParamBlockPtr);
  677. { Can be left emtpy if not needed. }
  678. { Called when calculating is through. See manual for more info }
  679.     begin
  680.     end;
  681.  
  682.  
  683.  
  684.  
  685. {***********************************************************************************************}
  686.  
  687.     procedure CleanUp (pb: ExtModulesParamBlockPtr);
  688.     { called when the function or program is removed from pro Fit's menus }
  689.     { in most cases, this function can be empty }
  690.     begin
  691.     end;
  692.  
  693.  
  694.  
  695.  
  696.  
  697.  
  698.  
  699. {***********************************************************************************************}
  700.  
  701. { This is the main procedure through which all calls to the external module go.                    }
  702. { Main takes care of calling the right procedure with the right parameters depending on        }
  703. { the value of "selector".                                                                            }
  704. { You don't need to touch this procedure                                                            }
  705.     procedure mainMain (selector: integer; pb: ExtModulesParamBlockPtr);
  706.     begin
  707.         Startup(pb);
  708.         case selector of
  709.             kSetup: 
  710.                 begin
  711.                     pb^.requiredGlobals := 0;
  712.                     pb^.versionNumber := VERSIONNUMBER;
  713.                     if sizeof(extended) = 10 then
  714.                         pb^.codeType := CPU68noFPU
  715.                     else if sizeof(extended) = 12 then
  716.                         pb^.codeType := CPU68FPU
  717.                     else
  718.                         pb^.codeType := CPUPowerPC;
  719.  
  720.                     SetUp(pb^.moduleKind, pb^.name, pb^.requiredGlobals, pb);
  721.                 end;
  722.             funcInitialize: 
  723.                 begin
  724.                     pb^.hasDerivatives := false;
  725.                     InitializeFunc(pb^.hasDerivatives, pb^.descr1, pb^.descr2, pb^.numberOfParams, pb^.a0, pb);
  726.                 end;
  727.             funcCheck: 
  728.                 pb^.answer := ord(Check(pb^.paramNo, pb^.a0, pb));
  729.             funcFirst: 
  730.                 First(pb^.a^, pb);
  731.             funcFunc: 
  732.                 Func(pb^.x^, pb^.a^, pb^.y^, pb);
  733.             funcDerivatives: 
  734.                 Derivatives(pb^.x^, pb^.a^, pb^.dyda^, pb);
  735.             funcLast: 
  736.                 Last(pb);
  737.             kcleanup: 
  738.                 CleanUp(pb);
  739.             otherwise
  740.         end;
  741.     end;
  742.  
  743. end.